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1.-  Introduction. 

Dolnf^  maximum  likelihood  estimation  when  the  like- 
lihood function  is  so  complicated  that  it  becomes  very  diffi- 
cult to  deal  with  it  ajialytically  ia  a  problem  which  has  to 
be  solved  numerically. 

The  process  of  estimating  the  parameters  of  an  un- 
derlying lognormal  distribution  when  sampling  is  considered 
to  be  without  replacement  and  proportional  to  size  constitu- 
tes one  such  case.  This  particular  sampling  process  is  being 
used  as  a  model  for  inferring  the  distribution  of  yet  undis- 
covered oil  and  gas  deposits  from  a  sample  of  discovered 
pools.  Pool  discovery  (observing  deposit  sizes  in  order  of 
discovery)  "is  more  akin  to  sampling  without  replacement  and 
proportional  to  size  than  to  sampling  values  of  independent, 
identically  distributed  random  variables"  (1).  This  feature 
had  been  ignored  until  recently. 

In  addition,  as  it  is  also  pointed  out  in  (1),  the 
population  from  which  the  discovery  process  picks  pool  sizes 
can  be  considered  as  a  sequence  of  N  values  of  N  mutually  in- 
dependent rajidom  variables  identically  distributed  with  a 
common  lognormal  density,  generated  by  nature.  There  is  em- 
pirical evidence  showing  that  a  variety  of  actual  size  dis- 
tributions can  be  reasonably  characterized  as  being  lognor- 
mal. 

The  mathematical  analysis  of  sucha  a  sampling  pro- 
cess is  contained  in  (1),  where  an  asymtotic  expansion  of 
the  likelihood  function  (valid  for  large  N  and  fixed  sample 
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size  n)  is  developed,  since  working  with  the  exact  expresions 
turns  out  to  be  very  difficult. 

In  turn,  since  this  asymtotic  expansion  is  analyti- 
cally complicated,  to  understajid  its  behavior  we  supplement 
the  treatment  of  it  in  (1)  with  numerical  analysis  of  parti- 
cular cases.  Its  complexity  may  be  better  understood  by  loo- 
king at  its  expresion.  The  expansion  for  the  sampling  densi- 
ty of  Y,  valid  for  large  p  =  N-n  isi 


^N 


n  r /  ^  ^T  X   n    Y.f  (Y .  [e) 


'J  {lt|pVg2(pMj^,Y)  +  ^P^V^g^(pMj^,Y) 


+  p(|M^-^^M2+-3-Mj)g^(pMj^,Y)  +  0  (p  ^)  } 
where  the  functions  g^(pMj^,Y)  are  given  by 

g^(pMwY)  =  I    (pM,+b.)"^  +  [  I     (pM  +b  )"^]^, 


1,3._  ^  ,  .,  .^  ,-l_  V'    „   ,  .-2. 


r  (pM  ,Y)  =  [  I     (pM,+b  )  ^]    +3[    I     (pM,+b  )  "-][    I     (pM.+b.)"'^] 

+  2  I     (pM  +b  )  ^, 

i=l      ^ 


+   3[    I     (pM  +b.)"l]^[  I     (pM,+b.)"2] 
j=l    ■"   ^       i=l    -^   3 


n  -  -,     n 

j=l    -      ^  j=i    1   3 


+  3[  I  (pM^+b  )  ^]^  +  6  j:  (pM^+b  )  ^ 


+  6[  I  (pM  +b  )  "■][    I    (pM,+b.)  -^l 
j=l    "•   =•      j=l    -^   3 


and  where  Y  is  a  vector  of  random  observations,  f(Y  |9)  is  the 

lognormal  density  with  vector  parameter  9=   (  /^  ,0"  ) ,  M^  and 

M^  are  the  first  two  moments  of  such  density,  N  =  population 

^  J 

size,  n  =  saunple  size  and  b.  =  ^  Y.  . 

Our  main  goal  in  studying  this  expansion  is  to  obtain 
approximate  maximum  likelihood  estimates  (MLE)  for  the  para- 

2 

meters  of  the  underlying  lognormal  density,  K  and  <J     ,    and 
for  the  finite  population  size  N.  Namely,  defining  the  like- 
lihood function  for  K  ,  CT^  and  N  given  a  data  vector  Y  as 
1(  A.  O-'^.nI  Y)  ,  we  wish  to  find  a  triple  (  K,,  (?lf  ,Nq)  of  va- 
lues (  k,  (7    ,N)  that  maximize  1. 

ML  estimators  may  be  analytically  obtained  in  the 
limit,  as  N  -»o<3  ;  hence  a  question  of  interest  is  how  large 
should  N  be  to  actually  obtain  in  practice  a  MLE  of  all  three 
parameters  J^  ,  Xr^   and  N.  It  is  also  interesting  to  get  an  idea 
of  how  stable  these  estimates  turn  out  to  be  as  a  function  of 
the  sample  size. 

The  approach  we  took  to  numerically  analyze  the  be- 
havior of  the  likelihood  function  for  particular  samples  was 
to  construct  iso  -  contour  graphs  of  it.  Doing  so  in  two  di- 
mensions with  one  parameter  fixed  allowed  us  to  visually  de- 
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tect  interesting  properties  of  the  function. 

What  we  did  was  the  following!  Taking  the  population 
size  N  as  fixed  and  regarding  \^   and  <r  as  arguments  of  the 
likelihood  function,  we  generated  sets  of  points  in  the  (  A,'^^) 
plane  satisfying 

1  (/^cr\  N  I  Y)  =  k, 
with  k  fixed,  (Such  sets  of  points  are  called  iso  -  contours). 
Plotting  sets  of  such  iso  -  contours  for  different  values  of 
the  population  size  and  different  values  of  the  sample  size, 
we  gained  some  insight  into  properties  of  the  likelihood  func- 
tion. 

In  what  follows  we  will  describe  these  properties  , 
along  with  the  estimates  we  obtained  for  the  parameters  men- 
tioned above. 

In  a  different  context,  certain  methodological  de- 
tails of  the  design  of  sin  iso  -  contour  generating  routine 
turn  out  to  be  interesting  as  well.  While  it  is  easy  to  find 
a  set  of  points  satisfying  an  iso  -  contour  equation  (over 
given  intervals  for  the  function's  arguments,  and  provided 
we  have  tabulated  function  values  for  points  belonging  to 
those  intervals),  what  constitutes  a  more  difficult  problem 
is  to  order  them  in  such  a  majiner  that  a  plotting  device  can 
go  from  one  point  to  the  next  actually  drawing  the  iso  - 
contour  curve.  Thus,  we  will  also  describe  the  method  we  used 
in  writing  the  routine  which  generates  ordered  sets  of  points 
belonging  to  a  given  contour.  Some  details  are  particularly 
relevant  when  one  is  interested  in  the  function's  extreme 
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points;  these  will  he  emphasized  in  the  corresponding  section. 
Program  details,  including  actual  program  listings,  are  inclu 
ded  in  Appendix  A, 

As  for  the  routine  implementation,  we  have  it  wor  - 
king  in  the  TROLL  system  environment,  taking  advantage  of  cer 
tain  plotting  facilities  already  available  in  this  system. 

A  couple  of  macros  which  make  the  routine  invocation 
easier  from  TROLL  are  included  in  Appendix  B,  along  with  their 
operating  procedures. 

Finally,  some  data  referring  to  program  requirements 
(memory  and  execution  times)  are  outlined  in  Appendix  C, 
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2.-  Iso  -  contour  p:eneratinf^  routine.  Outline  of  the  method 
employed . 

As  stated  in  the  previous  section,  the  main  diffi- 
culty lies  on  the  proper  ordering  of  the  set  of  points  belon- 
ging to  a  given  iso  -  contour.  The  procedure  we  followed,  ins 
pired  by  that  described  in  (2),  orders  the  sets  of  points  as 
it  generates  them  and  solves  other  problems  that  sometimes 
arise,  such  as  having  two  or  more  disjoint  iso  -  contour 
branches  in  the  intervals  of  interest  for  the  function's  ar- 
guments. 

Although  some  ])ublished  material  was  available  des- 
cribing programs  also  motivated  by  the  method  introduced  in 
(2)   ref.  (3),  (4)   ,  they  were  badly  documented  and  not  par- 
ticularly well  suited  to  interface  with  TROLL.  Further,  some 
interesting  problems  arise  when  generating  an  iso  -  contour 
near  the  function's  extreme  points,  which  were  precisely  the 
ones  we  v/ere  more  interested  in.  Our  final  procedure  differs 
slightly  from  that  in  (2)  regarding  to  these  problems,  amd  so 
we  thought  it  was  worthwhile  to  briefly  describe  it  here. 

A  general  outline  of  the  method  will  be  presented 
first,  so  that  those  problems  may  be  more  easily  understood. 

What  the  routine  accepts  as  input  data  is  a  table 
of  function  values  in  the  points  of  a  rectangular  grid  defi- 
ned by  a  set  of  rectangular  intervals  of  the  function's  ai^- 
guments,  with  increments  that  may  vary  along  the  intervals. 
In  such  a  grid,  two  units  turn  out  to  be  of  importance  when 
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the  aim  is  to  generate  iso  -  contours;  namely  the  grid  edges 
and  the  grid  cells. 

Grid  edges  are  used  as  units  during  the  first  step 
of  the  procedure,  in  which  each  one  is  checked  to  see  whether 
or  not  it  is  crossed  by  the  iso  -  contour  being  built.  A  way 
to  conduct  such  a  check  is  to  test  the  condition 


(f^  -  k)  .  (f2  -  k)  ;^  0 


(a) 


for  every  edge  in  the  grid,  where  f.  and  f^   are  the  function's 
values  at  the  edge  extreme  points  and  k  is  the  level  of  the 
desired  iso  -  contour.  Successful  tests  are  recorded  during 
this  first  step,  so  that  when  it  is  finished  (having  tested 
all  the  grid  edges)  ,  we  are  left  with  a  set  of  grid  edges 
crossed  by  the  iso  -  contour  under  construction. 

Two  properties  of  this  procedure  must  be  emphasi- 
zed at  this  point.  First,  test  (a)  will  only  detect  at  most 
one  crossing  point  at  any  edge  (  so  that  when  the  edges'  si- 
ze is  too  big,  situations  such  as  the  one  depicted  in  Fig. 
2.1  will  not  be  properly  recorded),  and  points  of  tangency 


Fig.  2.1 


between  grid  edges  and  iso  -  contour  lines  will  not  be  recor 
dcd  at  all.  Second,  the  procedure  will  go  on  in  its  second 
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step  to  locate  the  points  of  intersection  between  edges  amd 
contour  lines,  so  thkt  the  points  contained  in  the  resulting 
iso  -  contour  set  will  all  lie  on  grid  edges. 

The  second  step  of  the  procedure  takes  grid  cells 
as  working  units,  sind  relies  on  the  fact  that  when  the  grid 
size  is  properly  chosen,  either  none  or  two  edges  of  a  grid 
cell  will  be  crossed  by  an  iso  -  contour  line  (in  general,  if 
tangency  points  are  not  recorded,  an  even  number  of  edges  will 
be  crossed  in  any  grid  cell;  when  the  edges  size  is  too  big, 
again,  a  cell  may  be  crossed  in  its  four  edges,  this  is  a 
situation  likely  to  come  out  near  the  function's  extreme 
points  and  is  one  of  the  problems  that  will  be  analyzed  be- 
low) .  What  this  second  step  does  is  to  detect  an  iso  -  con- 
tour starting  point  (by  locating  any  crossed  edge  resulting 
from  step  1,  and  giving  preference  to  boundary  edges  to  pro- 
perly locate  possible  incomplete  iso  -  contour  lines  in  the 
region  defined  by  the  grid,  as  these  have  a  well  defined 
starting  point  at  the  boundary),  and  follow  the  contour  line 
through  adjacent  cells  to  the  one  containing  the  starting 
point.  This  operation  of  following  up  the  line  is  an  itera- 
tive one  which  is  done  as  followsi 

1.-  Having  located  the  edge  containing  the  starting  point, 
set  it  up  as  "current  edge". 

2,-  Delete  any  record  regarding  the  "current  edge"  as 
being  crossed  to  avoid  coming  back  to  it  thus  building  up 
endless  loops, 

3.-  Locate  an  iso  -  contour  point  in  the  "current  edge" 
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by  means  of  some  interpolation  procedure  between  edge  extreme 
points;  save  it  as  next  contour  point. 

4.-  If  a  "current  cell"  exists,  go  to  point  5;  otherwise 
decide  on  one  grid  cell  to  which  the  "current  edge"  belongs 
(in  general,  there  will  a  choice  among  two,  unless  that  edge 
is  at  the  boundary;  either  one  will  do) ;  call  it  the  "current 
cell". 

5.-  Of  the  four  cells  adjacent  to  the  "current  cell", 
pick  up  the  one  having  the  "current  edge"  in  common  with  it 
(notice  that  this  cell  will  have  either  none  or  one  crossed 
edge,  as  crossing  information  about  the  common  edge  was  pre- 
viously deleted);  call  it  the  new  "current  cell", 

6.-  Check  the  edges  in  the  "current  cell"  for  crossing 
information;  if  a  crossed  edge  is  detected,  call  it  the  new 
"current  edp;e"  and  go  to  point  2.  When  no  crossed  edge  is  de- 
tected, the  iso  -  contour  is  complete.  Record  this  fact  and 
go  to  point  1  to  look  for  other  possible  iso  -  contour  bran- 
ches. 

Such  an  iterative  process  will  be  better  understood 
by  following  the  sequence  of  sketches  in  Fig.  2.2,  where  (I) 
shows  the  actual  iso  -  contour  being  built  eind  successive 
ones  refer  to  procedure  steps  (a  x  stsinds  for  a  crossed  ed- 
ge, a  dot  for  an  iso  -  contour  point  in  the  resulting  set 
of  points;  the  shaded  cell  is  the  "current  cell",  the  darke- 
ned edge  the  "current  edge"). 

At  this  point,  we  are  in  a  better  position  to  dis- 
cuss the  kind  of  problems  that  are  likely  to  arise  in  the  neifh 
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Result  of  step  1 


I        M        I K 


(1) 


V7a 


y± 


(2,3.^) 


-T/ 


(5.6) 


(2,3.^) 


T/ 


(ZM^ 


(5,6) 


'ZIM- 


(2,3,4) 


is.e) 


7/ 


(2,3,4) 


Fig.    2.2 
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borhood  of  function's  extreme  points.  We  will  discuss  three 
of  them. 

As  has  been  emphasized  repeatedly,  grid  cell  size 
is  critical;  near  extreme  points  this  is  even  more  true. 

It  may  happen  that,  the  grid  size  being  too  large, 
an  small  iso  -  contour  near  a  function  extreme  point  will  be 
roughly  contained  in  one  grid  cell,  as  depicted  in  Fig.  2.3. 


Fig.  2.3 

Such  a  situation  violates  the  assumption  that  either 
none  or  at  most  two  edges  were  crossed  on  any  given  cell.  The- 
re is  no  way  of  deciding  in  which  way  the  four  points  in  the 
cell  should  be  joined  to  generate  the  iso  -  contour  curve.  The- 
re are  three  possible  ways  of  joining  them  as  shown  in  Figure 
2.4.  The  way  that  eventually  will  be  used  by  the  procedure 


Fig.  2.4 


described  above  depends  on  a  variety  of  factors,  such  as  on 
whether  it  picks  one  of  the  cell  four  crossed  edges  as  star- 
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ting  edge  or  not  and  the  implicit  sequence  on  which  cell  edges 
are  checked  for  crossing  information  upon  having  it  as  current 
cell. 

No  reasonable  way  oi  solving  this  problem  without 
changing  the  cell  size  (so  that  the  four  points  will  no  lon- 
ger belong  to  the  same  cell)  was  found,  unless  assumptions  are 
made  about  topological  properties  of  the  function  being  analy- 
zed. In  our  case,  it  is  not  safe  to  make  any  of  these  assump- 
tions -  at  least  in  principle.  In  addition,  given  that  one  of 
the  goals  of  our  study  is  to  determine  extreme  points,  a  finer 
grid  is  needed  in  any  case  to  obtain  precise  estimates.  Hence 
we  decided  to  apply  a  finer  grid  whenever  this  problem  appea- 
red. Curiously,  this  is  another  situation  which  points  out 
the  weaknes  of  software  procedures  in  dealing  with  pattern  re- 
cognition; most  of  the  early  approaches  to  such  problem  were 
very  similar  ;i.e.,  superimposing  a  grid  to  the  object  under 
inspection  and  being  unable  to  infer  deatails  from  the  gene- 
ral pattern.  In  our  case,  it  is  obvious  that  a  human  being, 
with  the  help  of  nearby  iso  -  contours  would  have  no  problem 
in  joining  the  four  points  in  that  cell  in  the  correct  way. 
The  procedure  doesn't  generate  such  information,  as  it  cons- 
tructs one  iso  -  contour  at  a  time,  but  there  is  no  apparent 

(iii 
way  to  make  "Mailable  to  the  routine  to  solve  the  problem 

unless  very  inefficient  prQoedures,  such  as  slope  matching 
are  used . 

A  second  problem,  particularly  relevant  when  the 
function  under  study  is  very  flat  near  the  extreme  points 
(and  this  was  our  case,  as  it  will  be  seen  later),  is  to  get 
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a  set  of  adjacent  grid  points  all  with  the  same  function  value, 
in  turn  equal  to  the  level  of  the  iso  -  contour  being  construc- 
ted. The  situation  is  illustrated  in  Figure  2,5.  Assume  the 
values  besides  grid  points  are  the  function  values  at  them, 
and  that  we  are  building  the  iso  -  contour  of  level  5.  What 
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seems  logical  to  do  is  to  record  the  edges  with  the  maximum 
value  at  their  extreme  points  as  being  crossed,  since  they 
themselves  are  so  closed  to  the  actual  iso  -  contour.  Test 
(a)  above  would  indeed  mark  those  edges  as  being  crossed,  but 
it  would  as  well  mark  all  the  other  edges  with  an  extreme 
point  of  level  5.  Then,  what  may  happen  (depending  upon  the 
same  kind  of  factors  outlined  for  the  previous  problem),  is 
that  we  may  end  up  with  an  iso  -  contour  plot  as  ridiculous 
as  the  one  depicted  in  Figure  2.6,  consisting  of  four  inde- 
pendent points  considered  as  iso  -  contour  disjoint  branches. 
Reference  (2)  proposed  a  way  of  avoiding  this  pro- 
blem by  disturbing  the  function  values  slightly  in  the  pro- 
per direction  so  that  only  the  outer  edges  would  be  recorded 
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as  crossed  (notice  that  when  those  points  are  at  a  relative 
minimum  the  disturbance  has  to  "be  in  the  opposite  direction)  . 


ir                                    . 

^^ 

Fip.    2.6 

1 

This  preserves  the  nature  of  test  (a),  but  search  for  such  grid 
nodes  has  to  be  undertaken  as  a  pre-step  for  the  method  in  use. 
The  way  we  solved  it  is  by  complementing  (a)  in  such  a  way  that 
edges  with  extreme  values  both  equal  to  the  iso  -  contour  le- 
vel are  not  recorded;  only  those  with  a  non  zero  difference 
at  one  of  their  extreme  points  are  so.  It  is  equivalent  to 
the  other  one  and  it  does  not  require  any  pre-step. 

Finally,  the  problem  depicted  in  Figure  2.1  is  one 
which  is  a  direct  consequence  of  the  cells  being  to  large.  It 
can  not  be  solved  by  einy  other  procedure  than  using  a  finer 
grid,  since  there  is  no  way  of  detecting  two  crossing  points 
in  an  edge  with  information  about  edge  extreme  points  (and  even 
if  it  were,  the  interpolation  process  would  end  up  with  the 
same  points  along  the  edges  generating  two  exactly  equal  bran- 
ches for  the  iso  -  contour. 
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3.-  Asymptotic  expansion  behavior  and  maximum  likelihood  esti- 
mation. 

We  had  at  our  disposal  two  sets  of  sample  data  to 
study  the  behavior  of  the  likelihood  function  generated  by 
themi  the  first  is  a  set  of  2k   North  Sea  pools  and  the  se- 
cond is  a  set  of  5^  pools  from  a  play  in  the  Western  Canadian 
Sedimentary  Basin. 

Before  going  any  further,  it  is  important  to  say 
that  the  chronological  order  in  reservoir  discoveries  is  of 
central  importance,  as  the  sampling  process  is  assumed  to  be 
without  replacement  and  proportional  to  -remaining-  size. 

Some  of  the  results  obtained  using  the  expression 
of  the  likeklihood  function  on  those  samples  are  summarized 
in  tables  1  and  2.  In  them,  the  values  of  ^   and   cr  which 
maximized  the  likelihood  function  for  different  values  of  N 
and  sample  size  n  are  shown  (When  the  sample  size  is  less  than 
the  actual  value  stated  above,  it  has  to  be  interpreted  that 
the  sequence  of  first  sample  values  was  used  as  data,  thus 
preserving  chronological  order) ,  The  estimates  for  N  =  ^^^ 
were  computed  analytically,  that  is, 

cr^=  —        ^       (log  Y,  -  log  g)^ 


A 


"    i=i 


log  g  -  ^  , where 

1 


,  /     log  Y. 

log  g  =  -_    ^1-      ^      J 


and  Y.  stands  for  the  j   observation. 

J 


TABLE    1 
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n 

N 

^o 

cr? 

Likelihood 
value 

Figure 

#_ 

2k 

200 

2.47 

1.^5 

-128.279 

1 

2k 

300 

2.45 

1.6 

-128.124 

2 

24 

303 

2.422 

1.569 

-128.12 

3 

2k 

400 

2.3 

1.7 

-128.179 

4 

2k 

500 

2.32 

1.6 

-128.269 

5 

2k 

600 

2.35 

1.6 

-128.369 

6 

2k 

750 

2.5 

1.3 

-128.463 

7 

2k 

999 

2.6 

1.2 

-128.314 

8 

2k 

2000 

2.6 

1.2 

-128.894 

9,10 

2k 

5000 

2.6 

1.2 

-127.551 

11,12 

2k 

10000 

2.6 

1.21 

-127.425 

13 

2k 

20000 

2.6 

1.21 

-127.362 

14,15 

2k 

io5 

2.59 

1.22 

-127.31 

16,17 

2k 

w) 

2.50286 

1.27988 

20 

300 

1.95 

1.95 

-103.827 

18 

20 

400 

1.75 

2.15 

-104.01 

19 

20 

06 

2.334 

1.31567 

12 

200 

0. 

5.2 

-   59.7028 

20 

12 

300 

0. 

5. 

-   60.0069 

21 

12 

oo 

2.1806 

1.48137 

TABLE    2 


WESTERN   CANADIAN   SEDIMENTARY   BASIN   DATA 
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n 

N 

h' 

crj 

Likelihood 
value 

Figure  # 

/ 

52 

500 

2.2 

2.3 

-314.024 

22 

52 

1000 

1.85 

2.5 

-310.084 

23 

52 

1500 

1.7 

2.6 

-308.58 

24 

52 

2000 

1.52 

2.7 

-307.813 

25 

52 

1.2648 

2.7213^ 

40 

1000 

1.85 

2.8 

-245.863 

26 

40 

1500 

1.5 

3. 

-244.921 

27 

40 

2000 

1.3 

3.2 

-244.442 

28 

40 

3000 

1. 

3.5 

-244.012 

29 

40 

1.40947 

2.95797 

30 

1000 

0.8 

3.9 

-183.915 

30 

30 

1500 

0.1 

4.7 

-183.3^2 

31 

30 

1.26892 

2.95797 

20 

1000 

-4. 

12. 

-126.765 

32 

20 

1500 

0. 

5.7 

-129.761 

33 

20 

2.02979 

2.74941 

3-^ 

The  column  labeled  "Figure  #"  indicates  the  iso  -  contour  plot 
corresponding  to  each  row.  In  these  plots,  the  axis  labeled 
"SIGMA"  actually  corresponds  to  cr  values. 

We  will  center  the  discussion  in  the  values  shown 
in  Table  1,    as  we  have  a  wider  range  for  them  and  the  operating 
characteristics  of  the  asymptotic  expansion  for  the  likelihood 
function  which  we  detected  seem  to  apply  to  the  values  on 
Table  2  as  well. 

The  main  feature  to  be  pointed  out  is  the  degree  of 
instability  which  shows  up  for  U    and   cr  estimates  as  N  stays 
constant  and  the  sample  size  is  chajiged.  More  concretely,  as 
the  sample  size  n  increases, ML  estimates  (keeping  N  fixed) 
for  U    and  t'' vary   systematically!  the  estimate  f or  JU.  in- 
creases while  that  for  cr  goes  dawn. 

Furthermore,  as  n  takes  values  below  a  certain  level, 
such  changes  in  yU  and  cr*"  estimates  become  more  apparent, 
taking  values  very  different  from  the  ones  obtained  for  lar- 
ger values  of  n.  This  is  the  case  for  n  =  12  in  Table  1  (even 
for  n  =  15  with  those  data  such  a  behavior  begins  to  be  very 
apparent)  and  for  n  =  30,  20  in  Table  2. 

There  are,  as  we  can  see,  two  two  main  reasons  for 
such  instability. 

The  first  is  probably  a  consequence  of  sample  da- 
ta variability,  which  is  relatively  more  important  when  the 
sample  size  is  small.  The  characteristics  of  the  samples  we 
had  -perhaps  due  to  the  fact  that  oil  producin;-:  areas  are  so 
recc'-nized  when  important  discoveries  take  place,  no  that 
some  early  observations  tend  to  be  large-  was  such  that 
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when  used  in  the  likelihood  function  produced  a  bimodal  beha- 
vior v/hich  was  hard  to  detect  as  the  second  relative  maximum 
builds  up  in  a  region  where  the  A  and  (T^  values  are  unreaso- 
nable as  estimates  of  the  underlying  density. 

Moroever,  the  second  relative  maximum  tended  to  be 
more  important  as  N  increased  for  fixed  n.  As  may  be  seen  in 
the  plots  corresponding  to  N  =  5OOO,  20000,  100000  when  n  =  24 
for  the  first  sample  (Figures  11,  12,  l4,  15,  16,  17),    it  rea- 
ches a  point  where  becomes  more  important  than  the  first  one. 
vVhen  this  happens,  those  values  out  of  the  reasonable  ranges 
for  A  and  <T     mentioned  above  begin  to  appear. 

For  very  small  sair!X)les  (e.g.,  12  in  Table  1,  20  or 
30  in  Table  2),  it  seems  that  the  second  relative  maximum  has 
taken  over  completely,  and  it  is"  the  only  one  obtained  in  tho- 
se cases. 

Before  discussing  the  other  reason  for  instability, 
m.ore  intrinsic  to  the  nature  of  the  problem  and  less  related 
to  data  variability,  it  is  perhaps  worthwhile  to  make  a  motho- 
dological  point  for  using  the  iso  -  contour  routine  or  a  si- 
milar tool  in  numerical  analyses  like  the  one  described  here. 
We  would  recomend  a  wide  range  on  the  function  arguments  to 
take  the  first  pictures,  forgetting  aboufreasonable"  values 
around  which  the  function  is  expected  to  behave  in  a  certain 
v/ay  (e.g.,  to  have  a  maximum).  Having  an  overall  idea  of  func- 
tion behavior  we  can  then  concentrate  on  certain  areas.  If 
this  step  is  omnitted,  the  first  results  tend  to  bias  the  ana- 
lysis to  pursue  a  certain  region  without  paying  any  attention 
to  others  which  may  contribute  heavily  to  explain  peculiar  be- 
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vior.  In  our  environment,  doing  so  implied  a  loss  of  detail, 
as  the  plots  we  generated  were  fixed  in  size,  but  even  in  this 
case  it  proved  to  be  in  the  right  direction. 

The  second  cause  of  estimate  instability  was  due  to 
the  extreme  flatness  of  the  function  near  its  maximum(s).  As 
it  may  be  seen  in  the  enclosed  plots,  this  circumstance  did 
show  up  rather  strongly.  Thus,  small  variations  in  the  func- 
tional form  due  for  example  to  small  changes  in  sample  data 
are  likely  to  produce  relatively  important  changes  in  the  lo- 
cation of  the  function's  extreme  point(s). 

The  changes  observed  in  the  location  of  maximums 
were  not,  however,  completely  arbitrary.  What  happened  was 
that  the  iso  -  contours  we  generated  were  roughly  elliptic, 
and  that  the  main  axes  of  such  ellipses  tended  to  conserve 
their  orientation  in  the  A  -  o^^  plane.  More  concretely,  the 
main  axis  tended  to  conserve  its  location,  so  that  different 
estimates  tended  to  lie  in  the  line  define  by  it.  This  may  be 
seen  as  a  common  characteristic  of  all  the  plots  we  include. 

As  a  function  of  N  with  n  fixed,  the  iso  -  contour 
plots  show  that  the  function  gets  tighter  along  the  small 
ellipses'  axis,  but  being  still  very  flat  along  the  other 
one.  This  suggested  that  the  underlying  colinearity  between 
U-   and   rr^  estimates  in  the  limit  as  N-*<»  (see  Jx.^  and  cT^ 
equations  above)  could  be  somewhat  preserved  for  finite  va- 
lues of  N.  Visual  inspection  of  the  plots  pointed  out  a  line 
equation  close  to 

h  =   3.775  -  0.8  cr^ 
for  the  first  sample  as  main  ellipses'  axis.  In  the  limit,  for 
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n  =   2k,    the  equation  Is 

^oo=   3.78274  -  crj  , 
not  very  far  off  the  one  visually  fitted.  This  would  explain 
the  overall  trend  of  JU  and  (T^  estimates  as  the  sample  size 
n  changes  which  we  described  at  the  begining  of  this  section  j 
it  may  help  to  foresee  their  behavior  for  different  sample  si- 
zes. 

We  were  also  interested  in  obtaining  an  estimate 
for  N,  the  finite  population  size.  Our  results  in  this  regard 
are  not  very  conclusive.  As  shown  in  Table  1,  consulting  the 
column  labelled  "Likelihood  values",  it  appeared  as  if  for 
N  =  303  we  obtained  an  absolute  maximum.  However,  as  we  had 
other  results  for  increasing  values  of  N,  the  function  value 
increased  above  the  quantity  -128.12  corresponding  to  N  =  303, 
and,  what  is  worse,  the  overall  maximum  (up  to  -116,511  for 
N  =  100,000  in  Figure  I6)  is  attained  at  a  point  with  unrea- 
sonable values  for  JU   and  cr^ .  Thus,  at  least  for  samples  as 
small  as  the  ones  we  had,  ML  estimators  for  N  appeared  to  be 
quite  instable,  and  other  estimation  methods  should  perhaps 
be  tried  to  ebtain  a  value  for  the  population  size.  It  has 
to  be  seen,  however,  how  the  asymptotic  expansion  works  with 
regard  to  N  for  samples  with  a  larger  n.  It  may  well  be  the 
case  that  our  sample  sizes  were  too  close  to  values  for  which 
the  second  function's  relative  maximum  begins  to  develop, 
thus  jeopardizing  the  function  usefulness  in  the  relevant 
range  for  ^-  and  cr^  . 

To  summarize,  the  two  reasons  for  estimates'  ins- 
tability discussed  above  have  pointed  out  how  sample  size 
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influences  the  estimates  we  obtain  with  the  analyzed  likelihood 
function  expansion.  Its  pathological  behavior  for  small  samples 
seems  to  suggest  that  having  more  numerous  samples  would  help 
to  obtain  better  estimates,  not  only  because  they  provide  mo- 
re information,  but  also  because  the  likelihood  function  begins 
to  behave  better  when  n  increases. 
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APPENDIX  A 
Iso-contour  plotting  routine  descriptlon_ and  progrram  1  i s t i ng^s. 

The  iso-contour  generating  routine  consists  of  a  set 
of  FORTRAN  and  370  Assembler  subroutines  which  form  a  package. 
CTOUR  is  the  name  of  the  main  routine  in  this  package,  and  its 
purpose  may  be  generally  described  as  follows i 

Given  a  tabulated  function  of  two  independent  varia- 
bles f(x,y),  it  determines  sets  of  points  belonging  to  the  ] o- 
cus  defined  by  the  equation  f(x,y)  =  C,  where  C  is  a  given  conj? 
tant.  (Subsequent  calls  to  CTOUR  with  different  C  values  allow 
to  obtain  a  set  of  such  loci). 

1 . -  Procedure . 

CTOUR  employs  as  main  procedure  the  one  described 
in  (l).  Roughly,  it  works  from  a  table  of  function  valuer;  at 
points  of  a  rectangular  grid  defined  by  two  auxiliar  input 
arrays  (which  permit  the  use  of  variable  increments  in  x  and 
y  along  their  respective  axes).  Given  a  value  for  C  in  f(x,y)=C, 
CTOUR  begins  by  marking  all  the  mesh  edges  crossed  by  the  lo- 
cus defined  by  the  above  equation  (an  auxiliar  bit  memory  is 
used  to  keep  track  of  edges);  then  it  computes  the  coordinates 
of  the  crossing  points  by  linear  interpolation  between  func  - 
tion  values  ad  adjacent  grid  points.  Points  belonging  to  the 
locus  are  returned  an  result,  ordered  in  Huch  a  way  that  it 
is  possible  to  actually  draw  the  contour  by  just  joining  thtMn 


(1)  Automatic  Contour  Map.  G.  Cottafava  and  G.  Lo  Roli,  Comm. 
of  the  ACM,  July,  19^9. 


A-2 

as  they  are  found  in  the  result  array.  This  result  array  may 
be  logically  divided  in  subtours,  when  the  locus  f(x,y)  =  C 
is  such  that  it  posseses  disjoint  branches  in  the  domain  de- 
fined by  the  grid, 

2,-  Input  and   output  variables  description. 

It  is  assumed  that  CTOUR  will  be  called  from  a 
FORTRAN  program.  Its  calling  sequence  is  as  follows i 

CALL  CTOUR ( A , NR , NC , XS , YS , ALEVL , lA , PTS , NDPTS , NPTS ) , 
where: 

A  is  a  bidimensional  array  containing  the  tabulated  func- 
tion values.  By  convention,  X  values  are  constant  along 
A's  columns,  while  Y  values  are  constant  along  A's  rows. 
For  convenience,  A  must  be  bordered  by  two  extra  rows  and 
two  extra  columns  which  do  not  need  to  contain  function 
values,  but  which  must  be  included  in  the  count  to  deter- 
mine A's  dimension. 
NR  is  the  number  of  rows  in  A,  as  declared  in  the  calling 
program  (including  the  extra  first  and  last  rows). 
NC  is  the  number  of  columns  in  A  as  declared  in  the  calling 

program,  also  including  the  2  extra  columns. 
XS  is  an  unidimensional  array  with  NC  elements  containing 

the  X  values  corresponding  to  different  columns  of  A. 
YS  is  another  unidimensional  array  with  dimension  NK  con- 
taining the  Y  values  corresponding  to  different  rows  of 
A,  (2) 
(2)  In  other  words,  A(I,J)  =  f  [XS(  J ) ,  YS(  I  )J  . 
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-ALEVL  is  the  function  value  for  which  the  contour  points  are 
to  be  determined. 
lA  is  an  auxiliary  memory  taking  the  form  of  an  unidimen- 
sional  integer  array.  It  must  contain  at  least 

((NR-1)*NC+(NC-1)*NR)/32|    elements,  and  it  must  be 
declared  accordingly  in  the  calling  program. 
PTS  is  a  bidimensional  array  where  the  locus'  points  are  to 
be  returned.  It  is  of  the  form  PTS(NDPTS, 2) ,  where  PTS(.,1) 
will  contain  X  coordinates  and  PTS(.,2)  Y  coordinates 
for  the  resulting  locus'  points. 

-NDPTS  is  the  maximum  number  of  points  that  may  be  allocated 

in  PTS  (i.e.,  its  first  dimension  as  declared  in  the  ca- 
lling program) .  If  CTOUH  reaches  this  maximum  during  the 
process  of  filling  PTS,  it  calls  another  subroutine 
(EXTEND)  which  is  supposed  to  save  the  contenlB  of  PTS 
when  enough  space  to  do  so  is  available,  cancelling  the 
whole  process  otherwise.  An  actual  program  for  EXTEND  is 
not  included  below  because  it  is  likely  to  vary  in  diffe- 
rent environments, 

-  NPTS  is  a  result  integer  variable  indicating  how  many  points 
are  returned  in  PTS. 

3.-  Example. 

In  what  follows  we  describe  how  to  set  up  the  input 
variables  for  a  call  to  CTOUR  in  a  concrete  case,  along  with 
the  output  p^enerated  by  it. 

Imagine  that  we  have  a  function  f(x,y)  tabulated  in 
the  interval   [l  6x^-3,  2^y^5]   ,  by  increments  of  0.5  in  both 
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X  and  y,  and  that  we  want  to  call  CTOUR  to  obtain  the  set(s) 
of  points  belonging  to  the  locus  of  level  4.5,  that  is,  be  - 
longing  to  the  curve  defined  by  f(x,y)  =  4.5. 

If  that  is  the  case,  the  input  data  should  be  orga- 
ni7,ed  as  depicted  in  Fig.  A.l,  and  the  grid  being  used  would 
be  the  one  shown  in  Fig.  A. 2. 

Assume  further  that  the  iso-contaar  f(x,y)  =4.5  has 
the  two  branches  drawn  in  Fig.  A. 2.  Then,  the  result  variables 
returned  by  CTOUR  upon  call  with  the  Fig,  A.l  variables  would 
be  of  the  form  depicted  in  Fig.  A. 3.  Notice  that  the  first 
point  of  a  contour  branch  is  repeated  as  its  last  one  to  ac- 
tually close  it,  unless  the  branch  is  incomplete  with  starting 
and  ending  points  at  the  grid  boundaries.  Different  branches 
are  separated  by  delimiters  labelled  NA  in  both  the  x  and  y 
coordinates. 

4.-  Brief  process  description,  with  references  to  actual  programs, 

All  subroutines  are  written  in  FORTRAN,  with  the  ex- 
ception of  MARK,  ERASE  and  SEE  which  manipulate  bits  in  the 
auxiliary  memory  lA  and  are  written  in  370  Assembler  for  impro- 
ved efficiency. 

b'ip-.    A. 4  shows  a  high  level  flowchart  for  CTOUH,  in- 
dicating calls  to  other  subroutines  in  the  package.  It  is  qui- 
te straight  forward,  although  some  steps  need  further  explana- 
tion. Following  are  succint  discussions  of  such  steps  and  the 
involved  subroutines. 

-  "Clear  auxiliary  memory"  just  sets  all  the  integer  ele- 
ments in  lA  to  zero. 

lA  is  used  as  a  bit  array,  each  bit  corresponding  to  one 
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Fig,  A.I.-  Input  variables  set  up  for  example, 
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Fip-.  A. 2.-  Grid  used  in  example.  f(x,y)='+.5  locus'  branches. 
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Fig.  A. 3.-  Result  variables  for  the  locus  in  Fig.  A. 2 

grid  edge.  Since  in  a  NC  x  NR  grid  there  are  (NC-1)*NR  +  (NR-1)»NC 
edges  and  one  370  word  holds  32  bits,  thus  the  assertion  about 
the  dimension  of  lA  in  point  2  above. 

Bits  in  lA  are  msLnipulated  as  follows  t  A  bit  is  set  to  1 
by  means  of  a  call  to  the  subroutine  MARK  whenever  its  corres- 
ponding grid  edge  is  crossed  by  the  locus  to  be  built.  Such  cir- 
cumstance is  checked,  for  every  edge  but  the  extra  ones,  by  the 
test 

(Al  -  ALEVL)  *  (A2  -  ALEVL)  ^0,         (3) 
where  Al  and  A2  are  the  function  values  at  the  edge  extreme 


(3)  Refer  to  section  2  for  test  details  in  special  cases. 
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Fig.    A. 4.-  CTOUR   flowchart. 
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points. 

The  correspondence  between  lA  bits  and  grid  edges  is  es- 
tablished as  follows  I  Edges  are  numbered  as  depicted  in  Fig.  A. 5 
for  the  case  of  a  4  x  7  grid  (vertical  edges  by  rows  first,  then 
horizontal  edges  by  rows).  Having  done  this,  given  an  edge  num- 
ber there  is  a  corresponding  bit  in  lA.  For  other  uses,  however, 
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Fig.  A. 5.-  Edge  numbers  in  a  4  x  7  grid. 

an  edge  is  more  conveniently  identified  as  shown  in  Fig.  A.6t  The 
vertical  edge  joining  the  grid  points  (I, J)  and  (I-H,J)  is  speci- 


(I.J) 


(I+1,J) 


(I,J+1) 


Fig.  A. 6.-  A  more  convenient  edge  specification. 


fied  by  means  of  the  3-tuple  (I,J,l),  and  the  horizontal  one 
joining  (I.J)  to  (I,J+1)  by  (I, J, 2).  The  subroutine  INDIA  trans- 
lates such  specifications  to  edge  numbers  to  manipulate  the  bit 
memory. 

-  The  subroutine  STPT  scajis  the  bit  indicators  in  lA  looking 
for  I's  by  mea^s  of  calls  to  the  subroutine  SEE,  which  accepts 

a  bit  number  and  returns  a  variable  =S  )  if  the  bit  is  <    >   . 
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Boundary  edpces   are  checked  first,  to  locate  incomplete  branches 
effectively.  STPT  returns  an  integer  variable  K  equal  to  the 
number  of  a  bit  in  lA  currently  set  to  1 .  If  no  such  a  bit  is 
found,  K  is  returned  set  to  zero. 

-  Calling  the  subroutine  EDGE  with  a  bit  number  K  returns 
an  edge  specification  in  the  form  already  seen  in  Fig.  A. 6,  by 
means  of  a  3-tuple  (I, J, KIN),  (i.e.,  EDGE  performs  a  transla- 
tion in  edge  specification  converse  to  that  one  done  by  INDIA). 

-  INTER  interpolates  linearly  in  the  edge  specified  by 
(I, J, KIN)  to  end  up  with  a  point  (XXX, YYY)  such  that 
f(XXX,YYY)  =  ALEVL  (with  the  error  proper  of  a  linear  interpo- 
lation; if  such  a  method  is  considered  ineffective  for  the  func- 
tion at  hand,  it  is  easy  to  write  an  alternative  INTER  subrou- 
tine to  perform  a  finer  interpolation,  although  this  circums- 
tance is  unlikely  to  arise  as  other  considerations  -see  section 
2-  force  the  grid  cells  to  be  of  a  size  where  interpolating  li- 
nearly is  usually  enough.) 

-  ERASE  accepts  a  bit  number  K  in  lA  and  sets  it  to  zero. 

-  NEXT  accepts  an  edge  specified  as  (I, J, KIN)  and  looks 
at  the  adjacent  cell  searching  for  another  edge  crossed  by  the 
iso-contour  under  construction.  If  such  an  edge  is  not  found, 
the  original  cell  itself  is  checked  for  a  crossed  edge  (becau- 
se on  the  boundary  it  may  happen  that  the  adjacent  cell  does 
not  have  another  crossed  edge).  When  this  search  also  fails, 
the  contour  branch  being  treated  is  complete. 

Cells  are  identified  by  means  of  their  north  -  west  cor- 
ner, and  their  edges  by  means  of  a  variable  set  to  1 ,  2,  3  or 
h   according  to  Fig.  A. 7.  Notice  that  the  k   edges  heve  to  be 
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identified  here  in  order  to  effectively  locate  an  adjacent  cell, 


(I.J) 


Fig,  A.?,-  Cell  and  edge  identification  for  NEXT, 

Given  a  cell  as  (I, J)  and  a  value  KIN(=1,  2,  3,  k)    indicating 
which  edge  was  crossed  in  that  cell,  NEXT  identifies  the  adja- 
cent cell  as  indicated  in  Fig.  A, 8  (where  x  indicates  the  cro- 

(I,J) 


KIN  =  1 


KIN  =  2 


(I 

J) 

— 

: 

■ 

— 

(I, J) 


KIN=3  KIN=4 

Fig.  A. 8.-  Adjacent  cell  identification  (NEXT). 


ssed  edge  in  cell  (I,J)  and  the  arrows  point  to  the  cell  con- 
sidered adjacent.  In  addition,  NEXT  updates  I,  J  and  KIN  to 
specify  the  next  edge  to  be  analyzed!  it  also  returns  in  K 
the  bit  number  in  lA  corresponding  to  it.  (K=0  indicates  no 
more  edges  crossed  by  current  locus'  branch,  thus  specifying 
its  end) . 
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■^'J     It-(  (  A(  1  ,  J  )-ALhVL)=M  A(  I  +  l,  J  )-ALFVI.)  )     6U,'^],70 
'5)      I  t^  (  A  (  I  ,  vJ  )-ALl-VL  )     6U,h^,60 
"■^      I  H  (  A(  I  +1  ,  J  )-ALhVL  )     (SO,  7  0,  6  0 

r  =;-  =;=  *  -;= 

r.=;:*.-;=*    V(-K"IICAL     hDGI-    CRdSSFD. 

C  «--=:= 

60     K  =  I  iMDl  A(  I  ,  J,  1  ,IMK  ,ImC  ) 

CAIJ.     MARK  (  I  A,K  ) 
70    Cum  1  INIie. 

r.:;=**.;=     (^KU-CT     A     CUMHMIK     SIARTUMf,     PlilMT. 

r  *  =:=  =;=  --^ 

7'^     r.  AIJ,     S  fP  )  (   1  A,MR  ,i\)C  ,K  ) 
IF(K)      no.  Ill), HO 

r,*=;:.;=:,:     AM     START  Il\)(,    RDIiMT     i-jAS     I- 1  llllMO- J  Oh'MT  1 1- Y     IM    WHICH    fcOt,^     IS     I 

c^.^=;=** 

t-'.O    CAIJ.     I-I)GH(K  ,I\|R,I\IC,  I  ,  J  ,K  f  M) 
C  -  *  =;=  ^;= 

r.  ;;:=;:*.;:     |)t-  f^-K|^/lJ|M^-     A|\)|)     SAVh     THH     C  HDR  O  I  i\J  A  I  h  S     l)F     THh     STARIJMI,     H-illM 
r. .;:;;:- ;;:     ^h      A  H  L  H      HI     CLilSh      II      WHt'M     bNIH-D, 

CALL     Il\ri  LR  (  I  ,.),K  Ii\i,XS,  YS,NR  ,NC  ,  A,ALFVL,  XXX  ,YYY  ) 
I-  [KS1  X  =  XXX 
I-  I  k  ST  Y  =  YYY 
(-.1  I      II  I     H6 


1  4  /  ;^  H  /  1  H 


C  =;=  =i'  =;=  >:<■ 


wo 


II)f-'l||FY     (,RIJSS|I\I(,     iMI)i\)|      IM     |-l)(,(-     t^Y     IN  I  l-.^Ml  11.  A  I  llllvl. 

CALL      I  IM  1  L  R  (   I  ,  J  ,  K  I  M  ,  X  S  ,  Y  S  ,  I\IR  ,  M  C  ,  /\  ,  A  L  h  V  L  ,  'C  X  X  ,  Y  Y  Y  ) 

FtiRi.FT     AHIIIII      IHIS    CRdSSlMG    PMI'Ml. 

CALL     FRASI-(  I  A,  K  ) 

IIPIlAT  F     Pilliv  I  S     SI-1  (  s  )  . 

i\l  H  I  S  =  w  H  I  S  +  I 

I  F  (  i-l)P  I  S-NR  I  S  )     H7  ,^-'0  ,40 
CALL     FX  I  l-i\iO  (  R  I  S  ,  iMOM  I  S  ) 
wk'!  s  =  l 

R  I  S  (  IM  R  I  S  ,  1  )  =  X  X  X 
R  I  S ( MR!  S ,/  )=YyY 


C  I  ll0i)4M() 
C  I  I  IO0S(u 
C  I  llfiKh  1  () 
CI  ll(W)hXn 
C  I  nf)Oh'-^() 
C  I  lUMi')A(' 
C  ]  IIOf)'^SO 
CI  I  lOOh^O 

c  I  iiO'isyo 
C1 noosHo 

C  (  (10(1  6  MO 
CT(U;060n 
CTf  lOOhlo 
CI  MOC'lS^O 
C  I  IIOil6-(0 

C  I  iifto^Ao 

(.  I  I  lOM^hO 
C  I  I  lOdhhO 
(.  I  I  lOOh  ^O 
C  I  IIO06Mfi 
C  I  (|006W() 
(.1  II0()  /()() 
C  1  110(1  /111 

c  inoii  /?o 

C  I  I  lOO  I  -iO 

(.  I  noo  /Ao 
C  I  I  io(i  ("•)() 

(,  I  IKMl  I  h[]  . 
(.  I  I  Kill  /  70 
C  I  I  lOll  /HO 
(.  I  I  lOl)  /Wli 
(.  I  (IO(iHOO 
(  I  I  lOOH  10 
C  I  II  O  0  M  ^  0 
C  I  IIOOM  -((I 
C  I  IIOOH^O 

C  iiinoHho 

C  I  I  lOOhM) 
C  I  I  KHiH  /I) 
C  I  I  lOOMHO 
(.  I  I  lOi  IMMI  I 
C  I  I  lOMHOd 

C  I  I  lOiiw  I  (1 
(.  I  I  i0(-w/O 

(.  I  I  KM  I'-J  -^O 
(,   I  I  KKIO'^.d 

C  I  I  lonwso 
(.  I  I  KMrvcid 


!l!l? 


A-14 


i-l  H/ASH     P.O 


CTIIDK 


D/iJ^^    =    7'Di20 


■:;==:^     IJKK     Kik    MhXI     CKIISS  Ii\i(,     MOINI, 
CALL     M  h-  X  7  (  I  /i ,  iM  K  ,  i\ii: ,  I  ,  J  ,  K  I  N  ,  K  ) 
f_ :;:.-;:  :;:;;:     | )  I  1 1- S      II      hXISi?      I  H     SO,       1 1  )  F  NT   I  H  Y      11      A  IM  I )     CiinTINIIF 

I  K(  K  )      ]  0(J,  lot),  H4 
r,v«:;=:;=     IF    iMilT,  THKMIimA  I  F    A     SIIHIHDK    L    \JU\K     IF    ANY     riTHFK     FXISIS, 

100  ( F( PTS( MPl S, 1 )-XS( ^ ) )      101,10^,101 

101  I  F  (  PTS(  NPl  S,  1  )-XS(  MCMl  )  )     10/>,10M,1U? 
1  o  ?  I  F  (  M  1  S  (  IM  P  T  S  ,  ;^  )  -  Y  S  (  /^  )  )     1  0  3  ,  1  0  w  ,  1  0  H 
lO:^  I  F  (  PI  S(  MP]  S,^  )-YS(  i'MRMl  )  )     lOA,  104,104 
lf)4  i\iPTS  =  NPtS+l 

IF ( WDPTS-MPTS )     104  ,106, 106 
106    CALL     FXTtWLM  PTS,imI)pTS  ) 

i\'  t-'  T  S  =  1 
106    PTSI MP  I S, 1 ) =F IKSTX 

PTS(  ivpTS,2  )=F1WS1  Y 
104    IMPI  S  =  iMPTS+l 

IF  (  rioPTS-MPTS  )     107  ,1UH,  1  OH 
1  o7    CALL     FXlFiMlX  p  I  s  ,i\iOPTS  ) 

i\i  P 1  S  =  1 
1  o  K    H  I  S  (  M  P  "I  S  ,  1  )  =  /  IM  /^ 
'~  P  TS  (  Ivp  rS,2  )=/iv)A 

(,i)      IM     76. 


TFl-iM  I  l\IA  I  F.      KtrlHKM      III     CALLI'Ml.     PKIICKAM. 


(•;  .■■  =;=  =; 

I 

TF 

I-;  i"i  I  i\i 

]  1 

o 

KF 
F:^ 

HIP  i\i 

|: 

1  6  / ;;'  3  /  I  M 

C  1  iioow/o 

( ,  1  1 1 0  n  w  H  ( ] 

C  1  1  KiOMMO 

(.  1  (I'M  ooo 

C  1 110  1  0  10 

(.1  iioio;;^o 

C  lll()l  0-(0 

C  1  iKll  O'tO 

C  1  lioio'^o 

CI  (101  0  60 

C  1(101  0/0 

C  1  110  1  OF,  U 

c  1  (10  1  04  0 

c  Mioi  yen) 

C  Kiol  1  10 

C  MiOl  1^0 

c  100  1  r-io 

C  1  (101  1  -^o 

c  r  ( 1 0 1 1  6  o 

(.  1  (  10  1    1   hi) 

C  1(1011  /o 

C  1  (K)I  1  H'l 

4 

C  1  (101  140 

C  1 1  (0  1  ;^oo 

(.  1  HO  I  /I  0 

C  1  (lol  ^/o 

( .  |(  1  0  1  /  H  o 

C  1  MO  1   XAO 

C  II  10  1  /'■.{) 

C  1  (10  1  /  U) 

(,  1  1  10  1  /  /o 

(.  1  1  (O  1  ^HO 

C  1  1  1  0  1  /■  4  0 

A-15 


:Lf:ASh  ?.i)  bUGE  UATH 

SUKKGDTINE  EDGE  (  K  ,i\)R  ;nC  ,  I  ,J,KIN) 
C   *-  TKAMSLATES  K->  I, J, KIN 

Q  *  if.  if  if 

KK  =  K-(  (  NK-1  )=;=|\IC  ) 

IF(KK)  1U,1U,?0 
10  K I M  = 1 

I  =  (  (  K-1  )  /|\iC  )  +  l 
J  =  K-(  (  I-l  )=;mMC  ) 
RETURN 
?0     I  =  (  ( KK-1  )/(NC-l  )  )  +  l 
J  =  KK-(  (  I-l  )*{KiC-l  )  ) 
KIN  =  ? 
RETURN 
END 


7b084 


10/31/07 


EOGOOOm 
El)(70()0?() 
E0(,000H0 
En(,00f)4() 

EO(,no()bo 

El)(,()0060 

EOGOooyo 

EDGOOOHO 
EOGOOO^O 
EOGOOIOO 
EOGOOl 10 
EDGOOl/^O 

Fn{;ooiHo 

EI)(,00140 
E0(,001  40 


t!*< 


I.t-ASF     ?.() 


INIHA 


DATh     =     7^083 


A-16 


11/3^/1  (S 


HI  INC  1  ION     I  Ml)  I 

:* 

COMHDTES     THb 

:;c  s;: 

K  Y     1  , J     L    KIN. 

♦  >l« 

IM     THE    CALL  IN 

11  =  1 

JJ=J 

KK  I  IM  =  K  IN 

(;()     1(1(40, 30»1 

10 

JJ=JJ+1 

K  K  I  N  =  1 

Gfl     Tn    40 

?() 

11=11+1 

K  K  I  M  =  2 

■^0 

INDI A=( NR-1 )* 

KPTORN 

40 

INOI A=JJ+( 1 I- 

Rl-I  IIKN 

fN\) 

A( I , J, KIN, MR ,NC ) 

IMDRX     IN     lA    CDRRH.SMdNIJING     Td    THE     hOGE     bPeCll-IFI) 

MR     t.    NC     ARF    NO.     OH     ROWS     i\    COLUMNS     OF     A     AS    DtCLARHI) 
G    PROGRAM! INCLUDING    MANDATORY     FXTRA    ROWS     AND    COLS.) 


0,  20)  , KK  IN 


NC+(NC-1 )-( I I-l )+JJ 
1  )^:=MC 


IiMDOOOl  0 
INI)000?() 
IMD00030 
INDOOOAO 
INOOOOSO 
iNDOOOhO 
INDOOOTO 
INDOOOHO 
INDOOO^O 
INDOOlOn 
INDOni 10 
IMDOOI 20 
IMDOO] 30 
INDOOl 40 
INDOOIHO 
INDOOl ^0 
INDOOl 70 
INDOOl HO 
INDOOIMO 
I  NI)002(ID 


n ! 

10  y 
I) 
p 
x 

x 

K 

10  I 

\) 

p 
y 

V 


A-iy 


[-ASF     ?.() 


INI  hk 


DAT'^     =     75079 


Ift/nO/A? 


SIIKKdDTiNK     INTI-K(  I,J,KIivl,XS.YS,NR,NC,A,ALFVL»XXX,YYY) 


10 


?() 

^0 


-^0 


CI  I 
k  I  I) 


NTIUJK     MDlNl      IS    (;bNEKATh(),V  I  A    LIMHAK     I  NT  bkPOL  AT  I(li\M  Y  I  NG     IN     THI- 
Ei)Oe     DbHINhO    HY     (I, J, KIN).     1  HE     kf^SMLT     IS     (XXX»YYY). 


II  MENS  I  ON     XS(MC  )  , YS (Nk  )  ,  A(Nk,NC  ) 
1  =  1 


(I     ni(  40,30,  10,?0  )  ,K  IN 

J=JJ+1 

I)     TO    40 

1=11+1 

YY=YS(  I  1  ) 

ELT  =  XS( JJ+I )-XS( JJ  ) 

=  ( All  I  , JJ  +  1  )-A(  I  1  ,JJ)  ) /DELI 

XX=(  ALEVI.-A(  I  I  ,  JJ)  )/k 

XX=XXX  +  XS( JJ  ) 

FTDkN 

XX=XS( JJ  ) 

ELT=YS(  I1  +  1)-YS(  [I  ) 

=  (  A(  I  I+l  , JJ  )-A(   1  I  ,JJ)  l/OELT 

YY=(ALEVI.-A(n,JJ))/P 

YY=YYY  +  YS(  I  I  ) 

ETDkN 

N|)    ■ 


INTOOOIO 
IMT00O?n 
I MT 00030 
iNl 00040 
IM I OOOSO 
INT00060 
IMl 00070 
INTOOOHO 
IMT  OOOMO 
IM 100100 
IMl 001 10 
INTOOl ?0 
IMT  00130 
IIVJTOOl  40 
IMTOOl SO 
IMl  001  ISO 
1  MI  001  /O 
[IM  I  001  HO 
1  Ml  00  I  MO 
I  M  I  OO/'OO 

1  M  I  (uvi  n 

IN  I  OO/'/^O 
INI  OO/'-JO 
I  MTOO>'40 
1  MVOO/'SA 


mf- 


m 

INIltX 

HUME 
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pa(;f 


SliDMCh-     STATHMEMt 


ASM     niO'D     ?1.1A 


1    h  K  fi  S  r 


STAKl 

FIMTKY  FRASI-^Sbe 

STM  1^,  1?,  i;'(l  3) 

RALK  KASF,0 

IISIiMf;  =;=,HA,SF 

F  0 1 )  K 

F  (M I        q 

Foil  10 

F  01 )  11 

F  (.)  U  ? 

Foil  3 

F  g  1 1        4 

I.  R(JMF,=F«1» 

LR  MASK.ROMF 

SLL  .     MASK, 31 
L  ADDR ,4(0,1  ) 

L  K ,0 (0, A DDK  ) 

LK  AC,K 

SK  ACkiinf 

SRL  ACS 

SLL  AC? 

LR  1NDFX,AC 

SLL  AC  3 

LK  6,  AC 

LR  ACK 

SK  AC'S  « 

SK  ACkiinf 

SKL  MASK ,0 ( AC  ) 

L  ADDR, 0(0, I ) 

L  ACO(  INDF  X  ,A|)DK  ) 

OR  ACMASK 

ST  AC()(  IMOFX  ,ADI)K  ) 

LM  14,  1?,  1? (  I  3  ) 

HR  14 

S  [  i-i  1  4  ,  1  ?  ,  1  /-  (  I   -i  ) 

HAIR        HASH,0 

IIS1N(t     «,hASF 

L  RIIMF  ,  =F  •  1  ' 

LR  MASK,  RIIMF 

SLL  MASK, 31 

L  ADDR ,4(0,1 ) 

L  K ,0( O, ADDR ) 

AC  ,K 

AC  ,R()MF 

AC  ,  S 

AC? 

INDFX,  Ai: 

AC,  3 

^,  AC 

ACK 

A  C  ,  /s 

A(,  ,RI 

1^1  A  S  K 


>:«GFT    K. 


.GFT     CDRRFSPdKiDlMG     ARRAY     ^-|FMF^M 


* 

«(,FT     HIT     POSITiniM 


.PRFRARF  MASK 

* 

*GFT  AFFFCTFI)  ARRAY  WORD 

MARK  SPFCIFIED  HIT 

STORF  WORO  HACK 
«RF  l.iIRN 


LK 

SR 

SRL 

SLL 

LK 

SLL 

LR 

LR 

SR 

SR 

SKI. 

1. 


IMF 

,0('\C) 
ADDR ,0(0,1) 
AC  ,D  (  1  iMDFX  ,  ADDR  ) 


»GFT  K 


GFT  CIIRRFSRONDIMG  ARRAY  FI.FmFMT 


*GI-1  H  I  1  ROSl  1  ION 
.PRhMARF  MASK 


=  (7F  T  AFFFC  1  FD  ARRAY  WDRD 


;iill«Ct 


A-19 


PAM- 


SliDKCb     STAll-MbMl 


ASM    Din's     ?1  .  lA 


XR  AC, MASK 

ST  AC  ,()(  INOFX  ,  ADDR  ) 

KK.  1^ 

S  T  M  1  <+  ,  1  ^  ,  1  ?  (  1  3  ) 

RALk  HASt-,{) 

l)SIN(-    «»KASh 

L  K()Mh,=F'l' 

LR  MASK,KI)lMl- 

SLL  MASK, 31 

L  Af)nK,4(0,  1  ) 

L  K»0  ((),  AODK  ) 

LR  AC,K 

SR  ACKilMH 

SRL  AC, 4 

SLL  AC,/ 

L  R  I  N I M-  X  ,  A  C 

SLL  AC,  -i 

LK  ft,  AC 

LR  AC,K 

SR  AC,  6 

SR  ACRIINf: 

SRL  MASK,0(AC) 

L  AODK, 0(1),  I  ) 

L  !?,()(  INDHX  ,A|-)I)R  ) 

L  A()I)K,H  (I),  I  ) 

MR  MASK,b 

BC  4,Yf-S 

SI  RdlMf- ,()(  0, /i|)|)K  ) 

R  HOT 

SLL  R(IMh,l 

SI  RnMI-,U  (0,AI)|)R  ) 

L  M  14,1^,12(13) 

KR  lA 

f^(\ll) 

=  1-  '  ]  ' 


HRASf^     SPeCIFIFI)    HIT 
STORI-    WflRI)    BACK 
*ReTl)RN 


«GFT     K 


(,f-T     CDWReSPONOlMC     ARRAY     hOMI-Mr 


«(;i- r  Bii    POSITION 

,HRf-HAR|-  MASK 

*(;t-f  AFFPCTHf)  ARRAY  WORO- 

« 

OT  NOYFS  ADDRFSS 

.CHECK  WHFTHFR  HIT  WAS  ON  OR  NOI 


*IF  NOT, SET  NOYFS 

«ANI)  RFTOKN 

.•IF  YFStSFT  NOYFS 

*RFTOKN 


TO,  1 
TO  ? 


<;iiH« 
ijIlK 

!l=l 

mfii 
ai' 

!FI» 

•«  JJ=. 

r,n  1 

fill' 

fill  ■ 
'J  11= 

UHI 
Ifll 

Kit 


KIM 


A-20 


MSF  ?.() 


MFXl 


DATf- 


7b079 


1S/4A/] 1 


SOBKODTINP  NtXi  (  lA,  MR, IMC,  I, J, KIM, K) 


10 

?n 

3(1 

AO 

40 

^0 
70 


LIHIKS 
(  I  ,  J  , 


90 


F(IK 
KIM  ) 


iMhXT     CROSSING    FDCtr,     Ih     AMY     C    kHUJKMS     IT     IN     HdTH    hf)RMS 
&    K.-KFhhrKR  IM(,     Td     I  A-/R  f:  TDK  MS     K=0     IF     IT     DdFSN'T     FXISI 


niMFNSIflN     IA(  1  ) 
I  1=1 
JJ  =  J 

FIIGF    DF     ADJACFNT     CFLL      ■ 

KKlN=KIN+2 

IF(KKIN-4)     2U,;^0,I() 

KKlN=KKIN-4 

(;f)    Tf)(  3U,A0,  bO,  60  )  ,KIN  i 

JJ=JJ-1 

GO     ril    70 

11=11-1 

Gil     III    70 

.Kl  =  .l.l+I 

(,ll     1(1    70 

11=11+1 

CAI.I.     CHFCKI  lA,  1  1  ,  JJ  ,KK  IM,MK  ,NC,K  ) 

1F(K)     80,80,91) 

WHFN  IT  DdFSN'l  FXIST  IM  '\  Ub     ADJACFMl  CFLL,CHFCK  THF  ClIkKhNI  (IMF 


CALL 
KF1IIK 
1  =  1  I 
,I  =  J,I 

KI^|=K 

RFTlIK 
F  M I ) 


CHFCK(IA,I,J,KIM,NR,MC,K) 
M 


K  IM 
M 


NFxoonin 

MFX000;^0 

MFxnoo3n 

MFXOOOAO 
MFXOnOhO 
MFXOOOAO 

MFxnooyo 

MFX00n80 
NF XOOOMO 

NFxnnioo 
MFxoni 1 0 

MFX001?0 

MFxnni 30 

MFxnOl AO 
MFXOOISO 
NFXOOl^n 

NF xoni 70 

MFXOOIHO 
NFXOOIMO 
MFXHO^-OO 
MFXOn?! 0 
MFXOn^^O 
MFXn(1?3() 
MFXOO^'AO 
MFXOO/'hO 
NFXOO^<SO 
MFXOO/^yO 
MFXOO^'Hn 
NFXno^WO 
NF  XOfV-iOO 
MFX0031  (1 
IMFXO0-i?O 
MFX0OH3O 

NFxoo34n 


If?. 


t?SCi 
If  BE' 
It 

1)1-1 

NCI 


|» 

Gil 

»IF 

1(1  w: 

.?  f,iii 

j=i 

IIH 
4LI 


A-21 


LFASF     ?.() 


SFARCH 


IJATR 


IbOHh 


0^/?^/n4 


U) 

?() 
?() 

40 


bO 


SKHKnilTIIV'F     SErAKCHi  A  ,MK  ,l\lCt^LFVL  ,(MCiMTK  ) 

SCANS     THfc     FABLF     A    AMI)    GhNEHATFS    MCNTK    VAl.UfcS    EtJDALLY     SPACfD 
HFTWf-FN     THF     MAXIMDi"!    AIM!)    MINIMUM    VALUBS     IIF    A. 

IIIMt-MSIfIN     A(  MR,MC  )  ,ALFVL(  NCN1K  ) 

AMAX  =  A(  1,1) 

A  M  I  i\i  =  A  M  A  X 

MRMl =MK-1 

MCM1=MC-1 

Dd    40     I=2tMKMl 

1)11     40     J  =  ?,NCM1 

IF( A(  I  ,  J)-AMAX  )     ^0,20,10 

AMAX=A( I t J ) 

(;(l     Til     40 

1  P  (  A(  I  ,  J  )-AMlM  )      -10,40,40 

A  M  I  M  =  A  (  I  ,  J  ) 

CIIMIIMOF 

j  =  Mf>jTK-l 

DFL  (  A=( AMAX-AM [ M  ) / J 

ALFVL  (  1  )  =AM  IIM 

DO     SO     1=1, J 

ALFVL(I+1)=ALFVL( I)+OhLIA 

KFl  llh(M 

FMI) 


Sf-AOOOlO 
SFAOOO^O 
SFAOnO-^0 
SFA00n40 
SFAOOO'jO 
SFAOOn^O 

SFAoonvo 

SFAOOOHO 
SFAOOOMO 
SF A  00 100 
SFAOOl  10 
SFA001?0 
SFA00130 
SFA00140 
SFAOOl bO 
SFA00160 
SFAOOl  fO 
SFAOOl HO 
SFAOOIMO 
SFA(M)2on 
SFAOO?! O 
SFAOO^/^O 
SFAOO?  -iO 
SFA00?40 
SFAOO^bO 


2.0 


STHT 


DATE  =  7b084 


14/00/17 


SUHKdUTINF  STPT (  IA,NK»NC,K ) 


LOOKS  FOK  AN  tl)Gh  MAKKtl)  IN  lA.IT  GlVt^S  PREFgRtNCE  Tl)  BUUNOARY 

ei)r7HS,sn  thai    incomplete   subtoors  are   located  eirst  and  their 

STARTING  POINTS  ARE  (JN  THE  HOONDARIES. 

DIMhNSION  I A( 1 ) 

NRMl=NR-l 

NCM] =NC-1 

1=?. 
b  Of]  10  J=2,NCM1 

K  =  INDIA(  I  ,  J,2»NR,NC) 

CALL  SEE( IA,K,NOYES ) 

GO  T(l(  10,90)  ,NOYES 
10  CONTINOE 

IE(  1-2  )  20,20,30 
20  I=NKM1 

GO  TO  b 
30  J  =  2 
40  DO  '->(.)     I=2,iMRMl 

K=INDIA(  I  , J,  1,NR,NC  ) 

CALL  SEE( IA,K,NOYES) 

GO  T0{ bO,yO ) ,NOYES 
bO  CONl INUE 

IE(J-2)  60,60,70 
60  J=NCM1 

GO    TO    AO 
70    NRM2=NRM1-1 

NCM2=NCM1-1 

DO    HO     I=2,NRM2 

DO    HO    J=2,NCM2 

f)0    HO    L  =  l,2 

K=IIMDI  A(  I  ,  J,L,NR,NC  )  -  • 

CALL     SEE( IA,K,NOYES) 

GO  r(l(  80,  VO  )  ,NOYhS 
HO  CONTINUE 

K=0  ■ 

90  RFTOKN 

END 
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STPOOOIO 
SI P00020 

srpooo^n 

STP00040 
SI POOObO 
STP0006n 

STpnno7o 

STPOOOHO 
STPOOOMO 
SI POOIOO 
ST  pool  10 
STP00120 
STP00130 
STP00140 
SI POOlbO 
ST  POO  160 
SI  POO  170 
STPOOIHO 
SI  POO  140 
SI  P0020() 
ST  POO?  10 
SI P00220 

SI  P002:-i0 

S  I  P0024() 
SI  P002b() 
S I P00260 
SI  PO0270 
SI P002H0 
SI  P0()290 
ST POOIOO 
ST  POO -5 10 
ST  POO:520 
STPOn330 
SIP00340 
STP00  3b0 
ST  P 003 60 
SI P003  /O 
STP()n3H0 
SI  P003M() 
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APPENDIX   B 

Routine  invocation  from  TROLL.  Macros. 

The  two  TROLL  macros  listed  below  facilitate  invoca- 
tions to  CTOUR  from  the  TROLL  system  environment.  (See  TROLL 

* 
manuals  for  details  in  specific  commands  ).  To  invoke  these 

macros  from  any  user  account,  it  is  necessary  to  specify  the 

following  command  after  logging  ini 

SEARCH  MIT128_MACR0| 

The  macro  e^tLIKPLOTS  includes  computation  of  likeli- 
hood function  values  at  grid  points,  so  that  it  is  only  useful 
for  our  particular  problem,  &PLOTS  is  more  general,  but  it  re- 
quires the  function  table  to  be  available  prior  to  its  invoca- 
tion. 

Calling  sequences  examples  are  also  included  and  are 
self-explanatory.  (Capital  letters  correspond  to  system  prompts) 


Available  from  the  National  Bureau  of  Economic  Research,  5^5 
Technology  Sq . ,  Cambridge. 


I.  IK^^|.||)S    - 
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KFMMdi^  r,(;(iTii   PKk- 

SHAKTH  HI  LL_FliiMr.  1  I  I  iM  ;  SFAKCH    w  IL'<_Fi)NCT  JlK-'  : 

llllTllPT  DFVTrp     T^^DIO; 

I  1 1  IT  H  P  T  M I  ]  I>/1  A  k  K  s  : 

(iiiTriMi    ^(i"»i;Riot    fiK    '  iMfi(,w  I  n  • :  "  : 

r.RFAIl     >,^v'iyiMIR     MM     V/AIJIFS:"     KHMD 

f.kFAn    (-.hoYiMiR    si-i;ma    so.    wfli.UFs:"    f,,Fi\m 

nil     KK4  =  FXPAwn(  K4,  IMDR  (  f.  5  )  »  KH  )H  (  f>  h  )   )  :  l)(i     K  f/i  =  I  k  A  l^' S  H  (  iO>  >  )  ; 

nil    K  s;  h  =  F  X  p  AMI  M  K  '^  .  ^!i  m  (  h  ^  )  ,  N 1 1 H  (  ^  A  )  )  : 

KRFAM     H'f'iH     VAt.DFS:"     f.FWD 

r>KFAn     ?;8"IIKSFRVAT  IIIMS  :  "     KFMO 

>'.KFA(1     >.4"M     VAI.IIF:"     NFMO 

0(1    f-,  «>;4  =  |.  IK!  HP  '  F  (  f.7,Kt.4,K't. '^  tCiiM-U  i\iF  (  KQ  )  ,  ^H  )  : 

DM      I  I- AK(.  (  1  )  =MAXS(  KHhy  )  ; 

KPRl'^'IA 

yilllK     MAXIMUM     VAI.IIF  IS     f.  I  F  AK(,-(   1   )   , 

r,Hll(iSI-     Yflliw     rOfvimiK     IFVFLS     AC,('.IIKI)TM(.LY  . 

K  hM  n 

?,  kFAI)     Kpi'MFW     LF\/FI_S?"      hhlMl) 

K 1 1-    r,;^    r,F(.)    i^wi    m,ii1(i    iild    f>iFFwi) 

DhlFTF     MAI  A     hX'iMFW     I.FVFI  S     nAi^if:": 

KSl  (.A: 

n (- 1) rr    K '■>.,],]  : 

ADD      I  1)1' ,  (•>()"FKI  I  hK      LFVFI.S      Ihl     MllllIrS:": 

F  I  I.F  : 

KMiT"    iH^-XT 
Kiil.n: 

Ki-iFAD      i,'-<"l|IJ)     l.h-\/i-I.S     KAMF:"      fihl^lD 
KhiFXT: 

nil   K  -^-Ki-y  (  siiR  1  ( ^  A,  f.'^ )  ) ; 

Dii    SI-  M  FAK(-(  f. '-^ ,  h-D  )  ; 

nil      KKK^r  =  MFDr,  I  1^  •  F  (  f.4,  (,  4,  ^H^<-)  ,  f,  -(  )    ; 

1)11    ?;•'■.  =^HK4r_x  : 
nil    K4.  =^Hf,'yr,_Y : 

1)11     I  F  A;vc,  (  ;^  )  =MliK(  ^,Hh4C  )  ;  , 

III  I    SF  M  F  AKr.  (  KHf, 9r. ,  4  )  ; 

^SFT  K  f  FAKf-,(  -^  )  =1  KFNI) 
KSF  T  K  I  h  ARM  A  )  =S  KFIMI) 
KSFTC.     KC,  I  F  AR(;(  I   )  =Ai-l(;i>FF(-HI  JKLl'II^H  IPOI-'S  I  IIVi-X  Y/      M- 1\<  I ) 

r, I.I  II  IDS  : 

f;SFT     K  I  F  AR(,(  1   )  =1      hFi\i|) 

n  F  I-  s  (^  A  r.  F    f,  A  .    K^  .  : 

KRFAD     ^;()'i\/ '\l'  I  Ahl    I-      IM     H(lR  I  /IIMl  AF     AXIS:"      -j-imD 
KTF      >;0     r.  FM      «Kh'>      ^(;l|■||)     RFVFS      KIH^-MD 

PFF'\iviF    1    ?: 

Kl^l  I  T'  I      D  I  Kl  I  I  X 

K  K  I-  \/  F  s  : 
ppi.  amf    ?    ]  : 

K,  D  I  H  D  I  X  : 

f,  S»-Tr     Kf  l.FAi-' (,(,>)  rrf.KHFPF     1      m;  I  FA.<;  (,(   I   )     M-ilD 

KSFir      >,r  t  FAI'M   1    )  =KS  I  R  I  PI        1      KC  1  I-  /u^i,  (   1    )      M-MD 

'Vi/\K  K      K  I  F  ARI-,  (   )    )      KC  I  F/\K(-,  (  ^^  )  : 

KSFT      K  I  F  ARM  1   )  ^>^  1  FAR(;(  h  I  FAi-;i;(  "  )    )+  I      M  iM  D 

K  IF      ^^  T  FAR  (;('-<  )      I.  I      K  I  FAKI-;(  ^  ) 

f.  S  F  T      f.  I  F  A  R  ';  (   -l  )  -  ^,  I  H  A  i--.  (.  (  ;-i  )  +  1       ):,  ^  I  I  D 

KSFV     >;  I  F  ARI- (  4  )  =r,  i  F  AR(;(  4  )  +  ]       f.Fl-lD 

yjw  n  II   D  1  Ki  1 1  X 

K  I  H  FiMD 

•sC,  '\  I  ; 

KPP  I  hi  I   ■\ 

/Mi.'.D /I-  (  !■■  S  I      -^     (   I  I'-'  llin-,      V'vl.ilFS      /"'I    (IMF     1^1-'.  I      A',      |i        I-  M    I        ,    A) 


»  )^   I     I 

»■:  I-  w  k  : 
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HI  J  ITS     - 
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M- K  K 1 1  k    K  ( ,  1 1  r  f  I    h  k  1^ 

SI    .''KTH     hIkST     i<  )  C  l-_Sf-(ii"ll-i^i  1  ; 
SH  /M<rH     Mill   _HI  liMC   I   I  I  Ihl  : 
1 1 M  1 1 1 M  T     I )  I-  V  I  ( :  (-     I  4  !  I  1  ( ) ; 
I  II  I  h  IP  t      iVl  n-i  /\kK  S  ; 

I II  i  I  •'  MM     K()'<  •  (,k  T II '     I  IK     '  i-M  i(,K  1  I)  •  :  "  : 
nil     1  H  Ak(-,(  1  )  =MAXS  (  M  "Yi  1111/     l/\Hl.i-:"): 
K  ^'  '•'  I  ivi  r  A 

I  HI-      ivi^X  (I'lllM     VAIJIH      I'M      i,\       IS      MH/\K(;(1). 
r".i-lilliSK     Yillik      rilMliilik     l.l-Vl-LS     AC(J  !'<;  I)  I  I  M.l.  Y  . 

Kh  ■■.Ml 

K'.'^AI)     ^,  ?'t|M(-l.-J      I,  P\/I-I.SV"      M-IMll 

K  I  h     K'^     f.hi-i     Kill     M,in  Ij     III.  n     M  hhh  I) 

1)1-1. f-[l-     DA  I  A     h'-^i'iMf-iAi     I.F\/I-I.s    MAi^il-:"; 

>,  s  T  (;  A  : 

1)1-1)11     f,  ^,  1  .  1  : 

'Will    Ti  Ik  .  K{i"i-iMl  hh-     l.l-vi-I.S     I  hi    oi  II I  I  HS  :  "  ; 

I-  M .  I-  : 

M-M|ii      iMf-Xl 

V,  I  I  I  ,  1  I  : 

MM- All    K'-<<iiiin    1. 1- \/ 1- IS    MAl-ii-:"    f,  hi-ili       ' 

K.  ■..  -  X  r  : 

1)1 .    h'^-kt-v  (  SI  IK  I  (  >,  -i ,  f,  -\  )  )  : 

III  I    Sf-"l  I  h  Aku  (  ^'-1 .  so  )  : 

M^i-;\n      >,4"Yi|llk      X      \//\|JI|-S:"      M-iMli 

>.i-''-Aii    V, h'l  Yi  !i  Ik    Y    v'Ai.uhs:"    M-i\iii 

Ki'i-A|)      >;(ii'Yllllk     r,l  li.'1  I  M  IK      r.lAMi-:"      M-i.||l 
I  M  I      f, 'S=:i-ihwr,  \i'  '  f  {  K'^,  Kh  ,  K]    .  t.  ■^)   : 

I II  ■    >>'■  .  =^h_y  ; 
III  I      K  '^  .  -  h(-i_  Y  : 

I II I     I  h  A'-'  (.  (  ?  )  =Mi  Ik,  (  Kf   )  ; 
I'll    sf-  M  I  'i)-;i.(  Ki'-,  '^  )  ; 

>. '-  I    I       K  i  I-  AK  I  .(-'.)-  I       M   lM|  1 
^  S'     I       >,  I  I    /Ik!  .  (  4  )  ='i      f,  I   ,Mtl 

>,si  I  r    Ml  1-  Ak(-,  (  I  )  --/\i,i,iii-i-(  li  I  JKl  .iivii  ikdi^s  1 1 1  \/ ivy  y  /    M-i-m 
(  I  M 1 1  n  s  : 

^M-  I        f,  1   h/\k'(,  (   1    1=1       ^hiMj  I 

i)H-i-  Sk  A. (M-     ^,4.     KS  .  ; 

?,i-;l-An      (-.il'iV '\k  I  AKI    I-       III     HI  IK'  I   /  I  11''^  A  I.      A  x   I  S  :  II      M-i,l|) 
K  n-      ^.f)     r.kO      "hh"      MJllll     Kl-VhS      ,',If--l\lli 
kMl.  Aim  I-      1       ■?  : 

;■;'.'  I  1 '  i    hi  Ki  n  X 

K I-'  k  \/  *-  s  : 

kki  A  ^'f-    >    ]  : 

M I  1  )<■  1 1  X  : 

f,  s  I-  I  r    f,  f  r  k  A  k  (M  >■ )  =  f.  K  k  k  jj  1     1    M ;  I  f-  /'.  <  ( .  (  1  )    k  k  !■  1 1 ' 

?,M-|C.      ?.f.  I  k'^-vf,!    1    )  -^  S  1  K  1  k|        1       ^(■.  I  k  'M-M    (    I)      M    I'll 

!■  A'-'K     >,  I  h  /\ki  .  (  1  )     M,  I  k/\Ki  ,(  ^  )  ; 

^Si    I       K  Ik  'il-M  .  (   1   )  =h  I  I    ni-l,  (  >,  I  k  /\k  (,(  'O   )+  1      M  I'll) 

^.  I  f.  1  I    AI.M  ,  (    M      I.   I       fM  »■  'll'l.(  /  ) 

■   ■-  I    1       .    I  I    /\k(,  (    /  )  =f,  I  I    Ai-M  ,  (    ^  )  +  I      f,  kl  'M 

>■'         1       t,  I  k  /MM  .  (  ''i  )  -  t.  1  I    AT  I.  (  'I  )  +  I       M-IMI  I 

M  .  I  I  M  I     I  I  I  k  I  !  1   K 

>,  I  I    I    -H' 

S  I    •\  I   ; 

^.  I '    ■    I  iM  I'  /I 

/  K  \  /       Ikik      kik.  S'l       ■<     (jii.iliMlk      v'ALllkS     Ai^kllHk      I-' k  N  I      /\  s      I  |. 

A:=>.  1  h-  A •-;!,(  h.i  )  ,  k  =  ^  1  k./\ -I,!  ',  1  )  .i,  =  h  I  k.'>)^i.(  ^,-'  )  . 
^  k-"l 
I  -Ml.!.: 


J, 1. 1  I    I  M       SI   (W\ 
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TROLL  rOf^MAND:  do  mu  =  seq(  0  .  ,  ^.  /  0. -L ) ; 
TROLL  rOMMANP:  do  s  i  ttma  =  seq(  0.  Oi,  3.  ,  0.  1 ) ; 
TROLL  rOMMANn:  Sdikplots 

'GRID*  OR  'fJOHRID*  :erid 

YOUR    MU    VALUtS:mu 

YOUR    SIGMA    SQ.     VALUf-S  :  s  i  ema 

B    VALUF-S:bb 

OBSF:RVATinNS:nsob 

U    VALUF;3U0 

-.YOUR    fiAXlMUM    VALUF-       IS    -lZi>.l25. 
niOOS^    YOUP    roMTOUR    L^Vf-LS    ArrORHINniY. 
fJf-U    LFVFLS?yps 

NF-.W    LFVf-LS    IIAMF:1pv300 

^.Nf-AV    SF-.RI  f-.S    LFV3UU 

hMTfcR    Lt-VF-LS    IN    OUOThS  :"-iZ8 .  xzij    -l^b.J    -IZb.o    -l^J    -UU    -ijj" 

VARIABLE    IN    ItORIZONTAL    AXISisima 


TROLL  rnrrANn:  Apiots  B-7 

'rRID*    OP     'NOHRin' :frrid 

YOUR    TAPLf- rnsobJUU 

T!^E•.    MAXIMUM    VALUt     IN    fJSORjUU    IS    -l^6.i^'j. 
CHOOSf-.    YOUR    CONTOUR    LfiVLLS    ACCOROINPLY. 
NtW    LhVl-LS?no 

OLP    LhVtLo    NAMf-:  leviUU 

YOUR    X    VALUtJ:mu 

YOUR    Y    V/"LUtS:slpma 

YOUR    roNTOUR    NAMh : nsob 5U0c 

'.VARIARLF    IN   HORIZONTAL   AXISrsirrma 
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APPENDIX  C 
Program  requirements  and  performance. 

The  memory  requirements  for  the  subroutines  inclu- 
ded in  the  iso  -  contour  generating  package  are  as  follows: 

Subroutine  Bytes 

CHECK  68k 

CTOUR  2,552 

EDGE  604 

INDIA  628 

INTER  1,134 

MARK,  ERASE,  SEE  240 

NEXT  822 

SEARCH  896 

STPT  1,088 


Total  8,648 

As  for  execution  times,  we  include  below  the  CPU 
times  taken  to  generate  and  plot  iso  -  contours  in  different 
cases.  We  measured  them  for  different  grid  sizes  (given  by 
the  number  of  grid  nodes)  sind  different  number  of  contours 
to  generate.  Actually,  these  times  are  execution  times  for 
the  macro  ("PLOTS,  and  thus  they  Include  some  TROLL  overhead. 
They  are  measured  in  seconds. 


C-2 
Grid  Size  Contours  frenerated   and    plotted 


25  X   25 

50   X  50 

15  X   15 

100   X   100 


3 

6 

0.67 

0.92 

1.42 

2.34 

2.43 

4.28 

3.82 

6.91 
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